﻿class correlador
{
private:
public:
	///Longitud de los correladores
	int ncorr;
	///Tiempo
	double *tiempo;
	///Desplazamiento cuadratico medio
	double *msd;
	///Correladores del calculo de la viscosidad
	double *A12,*A13,*A23;
	///Conteo para la normalizacion
	int *jcorr;
	///
	double *dn;
	///Tiempos de relajacion
	double *trelaj;
	///Almacenamiento posiciones particulas para msd
	double **x0c, **y0c,**z0c;
	///Valor para el tiempo de relajacion del primer canal
	double trelaj1;
	///Caja
	caja *box;
	///Paso de tiempo
	double dt, dth;
	///Determina la viscosidad aprox = (beta/6V)*lim t-> inf (A*A/t)
	double sumtau;
	///Solo componentes xz e yz del tensor
	double sumtau2;
	///Nombre del fichero
	parametros paramet;
	///Anchura de los campos de salida
	int ancho;
	///Almacenamiento de los tensores de esfuerzos para cada correlador
	double *tau0xy, *tau0xz, *tau0yz;

	correlador(caja *bx, parametros param)
	{
		ancho=20;
		paramet = param;
		ncorr=param.ncorr;
		box=bx;
		trelaj1=param.trelaj1;
		dt=param.Dt;
		dth=0.5*dt;
		///Contruimos los arrays
		tiempo = new double[ncorr];
		msd = new double[ncorr];
		A12 = new double[ncorr];
		A13 = new double[ncorr];
		A23 = new double[ncorr];
		jcorr = new int[ncorr];
		dn = new double[ncorr];
		trelaj = new double[ncorr];
		x0c = new double*[ncorr];
		y0c = new double*[ncorr];
		z0c = new double*[ncorr];
        tau0xy= new double[ncorr];
		tau0yz= new double[ncorr];
		tau0xz= new double[ncorr];
		for(int i=0;i<ncorr;i++)
		{
			x0c[i]=new double[box->npart];
			y0c[i]=new double[box->npart];
			z0c[i]=new double[box->npart];
		}



		inicializa();
	}

	~correlador()
	{
		///Eliminamos los arrays
		delete []tiempo;
		delete []msd;
		delete []A12;
		delete []A13;
		delete []A23;
		delete []jcorr;
		delete []dn;
		delete []trelaj;

		for(int i=0;i<box->npart;i++)
		{
			delete x0c[i];
			delete y0c[i];
			delete z0c[i];
		}
		delete []x0c;
		delete []y0c;
		delete []z0c;
		delete []tau0xy;
		delete []tau0yz;
		delete []tau0xz;
	}

	void inicializa(void);
	void calcCorrelador(double t);
	///Inicializacion del archivo de correlacion
	void initFileCorr(void);

};

void correlador::inicializa(void)
{
	for(int icorr=0; icorr<ncorr;icorr++) {
		dn[icorr]=1.0;
		jcorr[icorr]=0;
		A12[icorr]=0.0;
		A13[icorr]=0.0;
		A23[icorr]=0.0;
		if (icorr==0){
			///Valor del primer canal del tiempo de relajacion
			trelaj[icorr]=trelaj1;
		}else{
			///Se les da a todos los correladores un tiempo de relajacion enorme excepto al primero
			trelaj[icorr]=1.0e8;
		}
	}
}

void correlador::initFileCorr(void)
{
	ofstream MiArchivo (paramet.corrFile,ios::out);
	MiArchivo.precision(4);
	if (MiArchivo.is_open()){
		MiArchivo.setf(ios::fixed);
		MiArchivo<<"Correladores"<<endl;
		MiArchivo.fill('*');
		MiArchivo << setw(8*ancho) << ""<<endl;
		MiArchivo.fill(' ');
		MiArchivo<<"icorr"<<setw(ancho)<<"t-trelaj"<<setw(ancho)<<"dcmijx"<<setw(ancho)<<"dcmijy"<<setw(ancho)<<"dcmijz"
				<<setw(ancho)<<"sumtau"<<setw(ancho)<<"sumtau2"<<setw(ancho)<<"txy(0)*txy(t)"<<setw(ancho)<<"txz(0)*txz(t)"<<setw(ancho)<<"tyz(0)*tyz(t)"<<endl;
		MiArchivo.fill('-');
		MiArchivo << setw(8*ancho) << ""<<endl;
		MiArchivo.fill(' ');
		MiArchivo.close();
	}

}

void correlador::calcCorrelador(double t)
{

	double dd[3];
	double dcmij[3];
	double dcm;
	double dummy;
	int idnlast, idn;
    idn=-1;
	ofstream MiArchivo (paramet.corrFile,ios::out|ios::app);
	MiArchivo.precision(4);
    ofstream fsalida("posvel_bin.dat",  ios::out|ios::app | ios::binary);

	///Barre todos los correladores que ya han empezado a medir
	for(int icorr=0;icorr<ncorr;icorr++) {
		///dn(icorr) es el salto de tiempo del correlador icorr
		///Conforme el correlador tal va midiendo, lo hace con menos frecuencia cada vez.
		idnlast=idn;
		idn=int(dn[icorr]) ;
		///Si el tiempo por el que va el programa supera al tiempo de relajacion el
		///correlador entonces empieza a medir y se va contando las
		///veces que ese correlador aparece usando jcorr

		if(t>=trelaj[icorr]-dth) {
			jcorr[icorr]++;

			///Integrales 1,2; 1,3 y 2,3 por el metodo del rectengulo. alfa<beta
			///Se va integrando para calcular finalmente la viscosidad
			A12[icorr]+=(box->spxy+box->svxy)*dt;
			A13[icorr]+=(box->spxz+box->svxz)*dt;
			A23[icorr]+=(box->spyz+box->svyz)*dt;

			if (fabs(double(t-trelaj[icorr]))<dth) {
			    ///Primera medida del canal 'icorr'
				///Se asignan las posiciones del instante actual a las posiciones
				///de las que parte el correlador
				for(int k=0;k<(box->npart);k++) {
					x0c[icorr][k]=box->part[k].pos[0];
					y0c[icorr][k]=box->part[k].pos[1];
					z0c[icorr][k]=box->part[k].pos[2];
				}
				///Primera medida del tensor de esfuerzos del correlador icorr
				tau0xy[icorr]=(box->spxy+box->svxy);
				tau0xz[icorr]=(box->spxz+box->svxz);
				tau0yz[icorr]=(box->spyz+box->svyz);
			}

			///Si el numero de veces que ha contado el correlador es mayor que idn
			///escribe el tiempo por el que va menos el tiempo de relajacion del correlador
			///Se calcula el DCM para cada eje y se escriben las pos y vel de las particulas en
			///el archivo binario "posvel_bin.dat"

			if (jcorr[icorr]>=idn) 	{
//                A12[icorr]+=(box->spxy+box->svxy)*dt;
//                A13[icorr]+=(box->spxz+box->svxz)*dt;
//                A23[icorr]+=(box->spyz+box->svyz)*dt;
				double dt = t-trelaj[icorr];
				if(dt<0.) dt=0.;
				fsalida.write((char *)(&dt),sizeof(double));
				fsalida.write((char *)(&icorr),sizeof(int));
				dcm=0.;
				for(int k=0;k<3;k++) {
					dcmij[k]=0.0;
				}

				for(int j=0;j<box->npart;j++){
					///Escribe posiciones y velocidades
					for (int k =0;k<3;k++){
                        fsalida.write((char *)(&box->part[j].pos[k]),sizeof(double));
					}
					for (int k =0;k<3;k++){
                        fsalida.write((char *)(&box->part[j].vel[k]),sizeof(double));
					}
					#ifdef POLYDISPERSE
                        fsalida.write((char *)(&box->part[j].radio),sizeof(double));
					#endif

					///Correccion de la distancia entre las particulas i y j
					dd[0]=box->part[j].pos[0]-x0c[icorr][j];
					dd[0]-=(box->lbox[0])*anint(dd[0]/(box->lbox[0])); ///Entero + proximo manteniendo signo
					dd[1]=box->part[j].pos[1]-y0c[icorr][j];
					dd[1]-=(box->lbox[1])*anint(dd[1]/(box->lbox[1]));

					dd[2]=box->part[j].pos[2]-z0c[icorr][j];
					#ifndef CONFINED
                        ///Para la coordenada z no se aplican condiciones de contorno
                        ///periodicas cuando el sistema es confinado
                        dd[2]-=(box->lbox[2])*anint(dd[2]/(box->lbox[2]));
					#endif

					///Calculo de los desplazamientos cuadraticos medios
					for(int k=0;k<3;k++){
						dummy = dd[k]*dd[k];
						dcmij[k]+=dummy;
						dcm+=dummy;
					}
				}
				for(int k=0;k<3;k++){
					dcmij[k]/=double(box->npart);
				}
				dcm/=double(box->npart);

				///Determina la viscosidad aprox = (beta/6V)*lim t-> inf (A*A/t)
				dummy = (A13[icorr]*A13[icorr]+A23[icorr]*A23[icorr]);
				sumtau=(A12[icorr]*A12[icorr]+dummy)/3.;
				///Solo componentes xz e yz del tensor
				sumtau2=(dummy)/2.;
				///productos del tensor de esfuerzos correspondientes al canal y al instante
				dd[0]=(box->spxy+box->svxy)*tau0xy[icorr];
				dd[1]=(box->spxz+box->svxz)*tau0xz[icorr];
				dd[2]=(box->spyz+box->svyz)*tau0yz[icorr];

				if((t-trelaj[icorr])>1.0e-5){
					if (MiArchivo.is_open()){
						///escribe cuando ha pasado el tiempo de relajacion del correlador icorr
						MiArchivo<<icorr<<setw(ancho)<<t-trelaj[icorr];
						for(int k=0;k<3;k++) {
							MiArchivo<<setw(ancho)<<dcmij[k];
						}
						MiArchivo<<setw(ancho)<<sumtau/((t-trelaj[icorr]))<<setw(ancho)<<sumtau2/((t-trelaj[icorr]));
						//MiArchivo<<setw(ancho)<<sumtau<<setw(ancho)<<sumtau2;
						for(int k=0;k<3;k++) {
							MiArchivo<<setw(ancho)<<dd[k];
						}
						MiArchivo<<endl;
						//MiArchivo<<setw(ancho)<<sumtau<<setw(ancho)<<sumtau2<<endl;
					}
				}

				if (icorr<ncorr){
                    ///Si el desplazamiento cuadratico medio desde que se inicio el correlador es mayor
                    ///de 0.5 y el siguiente correlador no se ha iniciado, se inicia poniendo como
                    ///tiempo de relajacion el actual
					if ((dcm>0.5)&&(trelaj[icorr+1]==1e8)) 	{
						trelaj[icorr+1]=t;
					}
				}
				///Va incrementando el paso para escribir
				dn[icorr]*=1.25 ;
			}

			///Si se han cubierto todos los correladores y se ha superado un DCM de 25,
			///entonces se acaba la ejecucion
			if ((icorr==ncorr) && (dcm>25.)){
				break;
			}
		}
	}
	if (MiArchivo.is_open())
	{
		MiArchivo.close();
	}
	fsalida.close();
}
